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" ABSTRACT 

ro 

Aims. We present a method enabling the creation of constant-uncertainty/constant-significance light curves with the data of the Fermi- 
Large Area Telescope (LAT). The adaptive-binning method enables more information to be encapsulated within the light curve than 
with the fixed-binning method. Although primarily developed for blazar studies, it can be applied to any sources. 
Methods. This method allows the starting and ending times of each interval to be calculated in a simple and quick way during a first 
step. The reported mean flux and spectral index (assuming the spectrum is a power-law distribution) in the interval are calculated via 
" the standard LAT analysis during a second step. 

Results. The absence of major caveats associated with this method has been established by means of Monte-Carlo simulations. 
We present the performance of this method in determining duty cycles as well as power-density spectra relative to the traditional 
5^ ■ fixed-binning method. 

> ■ 
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i— i ■ 1. Introduction 

Variability studies of blazars provide a wealth of information on the dynamical processes at work in these objects and yield important 
I/"") . constraints on their physical parameters. Thanks to its all-sky coverage in 3 hours and large sensitivity, the Fermi-Large Area 
OO 1 Telescope (LAT), is enabling a long-term view of variability in the energy band 0.1-300 Ge V for a large sample of y -ray blazars. 
■"nJ" I About one thousand blazars are included in the Second LAT Active Galaxy Nuclei Catalog (lAcker mann et all 1201 lh . Due to the 
■ large variations in flux shown by blazars over different timescales, producing light curves that preserve the information provided by 
' the data is not an easy task. With a regular (fixed) time binning, using long bins will smooth out the fast variations assessable during 
fS| ■ bright flares. Conversely, using short bins might lead to hard-to-handle upper limits during low-activity periods. Here we present a 
t— I ■ simple method where the bin width is adjusted by requiring a constant relative flux uncertainty, cr\ n f = An, (alternatively a constant 
^ | significance) in each bin. This method allows for more information to be encapsulated within the light curve than can possibly be 
achieved for a fixed-binning light curve, without favoring any a priori arbitrary timescale and while avoiding upper limits. These 
advantages come at the expense of an iterative procedure required to find the proper bin widths satisfying the chosen criterion. The 
5— i \ elaborate, fairly computing-intensive maximum-likelihood procedure necessary to analyze LAT data does not lend itself very well 
to such an iterative procedure. While light curves with constant signal-to-noise ratios are commonly used at other wavelengths such 
as the X-rays, the above difficulty has so far hampered the generation of corresponding light curves in the LAT domain. The purpose 
of the method presented in this paper (referred to as the adaptive-binning method) is to provide a solut ion to this p roblem. 

The adaptive-binning method bears some resemblance to the Bayesian-Blocks (BB) method (IScargld H998). Both methods 
propose a time binning which adapts itself to the data instead of an a priori arbitrary regular binning. The main difference is that 
the BB method gives the most probable segmentation of the observation period into time intervals during which the photon arrival 
rate is perceptibly constant, i.e., has no statistically significant variations. A bin then defines the period over which a constant flux is 
observed. In the adaptive-binning method, the segmentation criterion is based on a given value of cr\ n F (or significance) regardless 
of whether the flux varies within the bin. The adaptive-binning method is expected to be more sensitive (i.e. allows the study of 
fainter sources) than the BB method since it makes use of both the photon spatial and energy information and properly accounts for 
the diffuse backgrounds while the BB method does not. The flare rise and decay times can also be determined more accurately with 
the adaptive-binning method. 

This paper presents the method and explores its possible inherent biases/caveats. For the sake of clarity, only light curves 
obtained with a constant relative flux uncertainty will be discussed, but changing the criterion to a constant significance, defined 
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from the Test Statistic^ (TS ), is straightforward. The method is described in Section 2. Validations with simulations are presented 
in Section 3. Section 4 presents results on timing analyses using the adaptive-binning light curves regarding the determination of 
the duty cycle and the power density spectra of the y-ray source. A summary is given in Section 5. 



2. Description of the method 

2.1. Optimization of the time bins 

The goal of the method is to compute light curves for an integral flux above a given energy E m j„, with constant Ao in all time bins. 
The source energy spectrum is assumed to be a power-law distribution with photon spectral index F. Although the method allows 
any E m - ln to be used, fluxes are reported in the following above the optimum energy (E{) for which the accumulation times needed 
to fulfill the above condition (i.e., bin widths) are the shortest relative to other choices of E m - m (see the Appendix for the derivation 
of Ei). The energy E\ depends on the signal/background ratio but is independent of exposure. It is calculated here using the flux 
and F determined over the whole time range spanned by the data. 

The method consists in solving the following equation for T\ (ending time of the interval) given a starting time To: 

<r blF (To,T u F,r) = A (1) 

where F and F are the average flux and photon spectral index over the interval [To, Ti] respectively. In order to avoid the penalty of 
excessive computing time required by the standard maximum-likelihood analysis, cr\ n F is estimated from the time-ordered list of 
photons contained in the region of interest (ROI), as described in the Appendix (see Eq. IA.8l ). 

Since the source is potentially variable, F must be evaluated consistently in each time bin. F is estimated through a procedure 
maximizing the likelihood (see Eq. IA.41 . this procedure being similar to a simplified version of the standard LAT analysis. In 
principle f should be evaluated similarly, but in practice f can be left fixedQ to speed up computation. Taking a constant F is largel y 
justified since spectral variations with time have been found to be very moderate in LAT-detected blazarfl (see lAbdo et a l. 2010b). 

In practice, due to the discrete nature of the data, Eq.[T]can only be solved in an approximate way. The procedure is the following: 

- For a given To, F (initially estimated over the whole period for the first bin) and f , the function cr ln F (To, Ti,F,T) is monoton- 
ically decreasing with increasing T\ (see Eq. IA.8b . The detection time of the earliest photon leading to the fulfillment of the 
condition cr ln F < Ao is hence taken as T\ . 

- F is then reestimated over [To, T\\ and <j\ n F is reevaluated using this new F value. Then 

• if (Ti n F is equal to Ao within a predefined tolerance, convergence is achieved and Tq is replaced by T\ ; 

• otherwise, T\ is reevaluated with the updated value of F. A bisection method, making use of the set of T\ estimates obtained 
in previous iterations, complements the procedure to speed up convergence. 

- The whole procedure is repeated until convergence is achieved. 

The times Ti for the different intervals are calculated sequentially until the condition <x ln F < Ao cannot be fulfilled using the 
remaining set of photons (these photons are left unused). Once the whole set of time intervals has been determined (this stage is 
referred to as "Step 1" in the following), F, T and cr\ n F are recalculated for each interval with the standard pylikelihood analysis 
("Step 2"). The consistency between these values of cr ln F and Ao is then checked. Using Monte-Carlo simulations, an excellent 
agreement has been observed in most cases, making further iterations unnecessary. 

Being driven by the number of accumulated source photons, the bin widths provided by this method depend on both the source 
flux and the exposure of the instrument for the source location. The modulation of the LAT daily exposure over the precession 
period (53.4 days) ranges between a few percent and +50% (for a rocking angle of 50°) depending on the source declination. 

2.2. Pylikelihood analysis 

Results described in the following were obtained via an unbinned pylikelihood analysis performed with the Fermi-LAT ScienceTools 
software packag^ version v9rl9p0. The tools used in Step 1 and the pylikelihood analysis made use of the same list of photons with 
energy above 100 MeV and contained within a 20°-in-diameter ROI. The P6_V3_DIFFUSE set of instrument response functions 
(IRFs) was employed. The Galactic diffuse emission model version "gll_iem_v02.fit" and the isotropic (sum of residual instrumental 
and extragalactic diffuse) background "isotropicjem_v02.txt" (files provided with ScienceTools) were used throughout the proce- 
dure. The normalization of both diffuse components were left free in the fit for the pylikelihood analysis. For uniformity over many 
realizations produced under the same conditions, the optimum energy was evaluated from the true flux and spectral index used in 
Monte-Carlo simulations. 



1 The Test Statistic is defined as TS = 2(ln X(source)- In _£(nosource)), where £. represents the likelihood of the data given the model with or 
without a source present at a given position in the sky. See Appendix for details. 

2 T can be taken from a source catalog when available (for real data) or set to the measured value over a long period. 

3 Preliminary light curves of bright LAT-detected blazars computed with real data proved that taking a constant photon spectral index in the bin 
width estimate provides reasonably accurate results in terms of <r ]n F . 

4 http://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/ 
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3. Validation with simulations 

In this section, we first present the results of the adaptive-binning method applied to simulated variable sources. Then, we address 
some of the possible caveats of this method by considering simulated constant sources. For this purpose, we compare results 
obtained for adaptive-binning and fixed-binning light curves regarding the following issues: biases and fluctuations in the measured 
fluxes and photon indices, correlation between flux and F and finally the flux correlation between consecutive bins. Fixed-binning 
analyses were carried out with the same total number of bins as the adaptive-binning analyses to enable meaningful comparisons. 
The criterion adopted for all adaptive-binning analyses is Ao=25%. 

3.1. Description of simulations 

Data were simulated in the energy range 0. 1-200 GeV with the tool gtobssim, part of the Fermi-LAT ScienceTools software package, 
which generates photon events from astrophysical sources and models their detections according to the IRFs. The P6_V3_DIFFUSE 
set of IRFs was used consistently with the analysis described above. These simulations do not include all known systematic effects 
of the real data, in particular related to the varying instrumental background during the orbit. However, this lack is not believed to 
signif icantly affect the conclusions drawn thereafter due to the moderate magnitude of these effects (a few % in the flux, lAbdo et aD 
2011). The simulations were performed for an 11-month long period starting at the beginning of the scientific operation of the 
Fermi mission, using the actual LAT pointing history. Although light curves shown in the following use fluxes above the optimum 
energy, the different examples are distinguished according to their true fluxes above 100 MeV (Fioo in 10~ 7 ph cm~ 2 s _1 ) to make 
comparisons clearer. All sources were simulated at the same high-Galactic latitude sky position, (l,b)^(80°, 80°). 

3.2. Light curves of variable sources 

In this section, the method is illustrated for three variable sources which were simulated by means of templates showing realistic 
flux variations typical of blazars. The source spectra were modeled with a single power-law function with an average flux Fioo=2 
(corresponding to a fairly bright blazar detected by the LAT) and a fixed photon index T-2.4. The optimum energy is £1=230 MeV 
for these parameters. 

FigureQ]compares the true light curves with light curves calculated using the adaptive-binning method for the three sources. The 
analysis for the light curve shown in Figure[T]top was also conducted assuming a reversed arrow of time, the results being presented 
in Figure [2] Statistically significant fluxes can be determined over very different time scales, ranging from a few hours to several 
months in these examples. The occurrence of significant short-term flaring activity during the low-flux states can be ruled out. 

An adverse feature of the method is the flare truncation (the light curve in FigureQ]: shows such a feature at MET^255000000). 
When a significant flare occurs after a long period of low activity, a non-negligible number of photons from the flare has to be "con- 
sumed" to counterbalance the accumulated background photons and enable the criterion to be fulfilled. Possible ways to mitigate 
this effect are being explored and will be implemented in future versions. The tools allow light curves to be generated by considering 
either increasing or decreasing photon times, so the robustness of any analysis results (time lags, duty cycles, rise and decay times...) 
can be tested by using both flavors. 

Figure [TJ; may convey the impression that the adaptive-binning method will lead to light curves lagging the real one. However, 
cross-correlation analyses have shown that no significant lag is observed for the light curves investigated here, probably because 
truncation affects only very few bins and the flare leading edges are still well sampled. 

Figure [3] shows the distribution of relative flux uncertainties calculated from the pylikelihood analysis for the three variable 
sources. All three distributions are centered at a value close to the target Aq—25% (depicted by the dashed line) with essentially all 
values included in the range 20%-30% with a typical rms of 3%, demonstrating the performance of the method. 

3.3. Caveats explored with steady sources 

Steady sources with three different fluxes (Fioo=0.5, 1 and 5) and two different photon spectral indices (r= 2.0 and 2.4) were 
simulated over the same 1 1 -month long period as for variable sources. Fifty realizations were performed for each case. The average 
numbers of bins in the light curves as well as the optimum energies are reported in Table[T]for the different cases. For Fioq=1 and 
T=2.0, a light curve where fluxes are estimated after Step 1 (i.e., using the simple maximum-likelihood procedure described in the 
Appendix) for a particular realization is compared to that resulting from Step 2 (i.e., where fluxes are computed with the whole 
pylikelihood analyses with the time intervals established in Step 1) in Figure |4] There is a good agreement between these fluxes, 
which holds for all realizations and parameter sets. Examples of adaptive binning light curves for the flux are shown in Figure|5]for 
the different cases, while light curves for the photon spectral index are displayed in Figure|6] FigureQcompares the flux distributions 
obtained via the adaptive- and fixed-binning methods by summing up all realizations. For a given set of conditions, the distributions 
have similar means and rms but different skews. A positive skew (typically equal to 1 cr for the distributions shown in Figure|7]) is 
a feature of the adaptive-binning method. A positive fluctuation in flux leads to a shorter-than-average bin associated with a high 
apparent flux. A negative fluctuation has less chance to be recorded since the bin will be extended until enough photons have been 
accumulated to meet the Ao-wise criterion. The photon index distributions have been found to be very similar for the two methods 
and close to gaussian distributions. 

The distributions of relative statistical uncertainties obtained with the adaptive-binning method (resulting from Step 2) are 
presented in Figure [8] These distributions are more gaussian-like than those obtained for variable sources, with a mean value close 
to 25% as required and a typical rms of about 1.7%. Simulations performed for a very low flux of F ioo=0.1 (r=2.0), where the 
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240 245 250 255 260 265 270 

MET (s) 

Fig. 1. Simulated light curves (blue) with the adaptive-binning light curves superimposed (Step 1: black open squares, Step 2: red 
solid circles) for three variable sources. MET stands for Mission Elapsed Time. 

method yields a single bin over the 1 1 -month period, show that the mean of the distribution remains within 10% of the target value 
even in that extreme case. 
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Fig. 2. Same as Fig. QJ but the analysis was conducted assuming a reversed arrow of time. 
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Fig. 3. Distribution of relative flux uncertainties given by the pylikelihood analysis for three variable sources. The dashed line depicts 
the target A =25%. 





r = 


2.0 


r = 


2.4 


^100 


< N bins > 


Ei (MeV) 


< N U ns > 


Ei (MeV) 


0.5 


14 


396 


8 


299 


1 


37 


332 


23 


253 


5 


268 


246 


192 


205 



Table 1. Average number of time bins and optimum energies for the different simulation conditions. 



The possible biases induced by this method are explored in the following. 



3.3.1. Bias and statistical fluctuations in flux and photon-index measurements 

The pylikelihood analysis is expected to provide the correct mean flux and photon spectral index over all possible time intervals. We 
nevertheless check for the absence of biases in the measurement of these parameters when the time intervals are determined with the 
present method. Furthermore, possible additional fluctuations in the measured parameters could arise from the intrinsic correlation 
between bin width and photon contents of the bin. Both effects are addressed here. 

One defines AF (AF) as the difference, expressed in sigma, between the measured flux (index) and the average value measured 
over the whole 11 -month period. Summing over all realizations, the distributions of the means (left) and rms (right) of the AF 
(top) and Ar (bottom) distributions are displayed in Figure|9]for a particular set of parameters. The results of the adaptive- (solid) 
and fixed- (dashed) binning methods are compared in this figure. The means < AF > and < AT > are listed in Tables [2] and [3] 
respectively, while the rms values are given in Tables|4]and|5] 

A small bias is observed in the mean flux values. Since similar features are seen for both adaptive- and fixed-binning light 
curves, one can safely conclude that this bias arises from small systematics effects pertaining to our pylikelihood analysis and is not 
created by the adaptive-binning approach. 
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Fig. 4. Comparison of flux from the tool, i.e., after Step 1 (open diamonds) and flux resulting from the pylikelihood analyses, i.e.. 
after Step2 (filled circles) for a simulated constant source with F ioo=l and F-2A. 
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Fig. 5. Examples of light curves for the different sets of parameters : F-2.0 (left) and T=2.4 (right), Fioo=5, 1, 0.5 from top to 
bottom. The solid line depicts the true flux value. 



^100 


Adaptive (T=2.0) 


Fixed (T=2.0) 


Adaptive (T=2.4) 


Fixed (T=2.4) 


5 


-0.20±0.01 


-0.15±0.01 


-0.19±0.01 


-0.15±0.01 


1 


-0.13±0.02 


-0.09±0.02 


-0.14+0.03 


-0.13±0.03 


0.5 


-0.06±0.03 


-0.07±0.03 


-0.12±0.05 


-0.13±0.05 



Table 2. Comparison of < AF > for adaptive and fixed binnings. 



The average rms values are found to be close to the expected value of 1 for all cases (with slightly lower values for the adaptive- 
binning approach relative to the fixed-binning one). One can thus exclude that the adaptive-binning approach gives rise to significant 
fluctuations in the measured parameters. 



3.3.2. Flux-index correlation 

Another issue worth investigating is the correlation between measured photon index and flux in a given time interval. All photons do 
not contribute equally to the fulfillment of the criterion, with higher-energy photons contributing more. The detection of a series of 
photons with higher-than-average energies will thus lead to the criterion being satisfied faster and consequently to the interval under 
consideration being closed earlier. In order to evaluate whether the adaptive-binning method creates additional correlation between 
the photon spectral index and flux (as measured in Step 2), the corresponding correlation factor is computed for each realization. 
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Fig. 6. Temporal evolution of the measured photon spectral index: F=2.0 (left) and F-2.4 (right), Fioo=5, 1, 0.5 from top to bottom. 
The solid line depicts the true photon index. 




Normalized Flux Normalized Flux 




Normalized Flux Normalized Flux 



50 : 




Fig. 7. Normalized flux distributions for adaptive (solid) and fixed (dashed) binnings and steady sources with r=2.0 (left) and T=2.4 
(right), Fioo=5, 1, 0.5 from top to bottom 



Since E m - m is chosen as the optimum energy, for which the correlation between photon index and integral flux is minimum by 
definition, the observed correlation factor is expected to be small. The correlation-factor distributions obtained for adaptive and 
fixed binnings are shown in Figure [10] For all cases, the distributions are very similar for the two approaches and are all centered 
at small values, with larger fluctuations for lower fluxes. Therefore, the measured photon spectral indices and fluxes do not appear 
more correlated when using an adaptive binning relative to a fixed binning. 
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Fig. 8. Distribution of cr\ n p for Arj=25% (depicted by a dashed line) for the different parameter sets: F=2.0 (left) and F=2.4 (right), 
^ioo = 5, 1, 0.5 from top to bottom. 



^100 


Adaptive (r=2.0) 


Fixed (r=2.0) 


Adaptive (T=2.4) 


Fixed (T=2.4) 


5 


0.02±0.01 


0.02±0.01 


0.02±0.01 


0.01+0.01 


1 


0.02±0.03 


0.05±0.03 


0.04±0.03 


-0.01+0.03 


0.5 


0.06+0.03 


0.05±0.03 


0.03+0.06 


0.03+0.05 



Table 3. Comparison of < Ar > for adaptive and fixed binnings. 



^100 


Adaptive (T=2.0) 


Fixed (r=2.0) 


Adaptive (T=2.4) 


Fixed (T=2.4) 


5 


0.99+0.01 


1.06+0.01 


0.97+0.01 


1.03+0.01 


1 


0.96+0.02 


1.01+0.01 


0.93+0.02 


0.98+0.02 


0.5 


0.90+0.03 


1.02+0.03 


0.81+0.03 


0.87+0.03 



Table 4. Comparison of the average rms of the AF distributions for adaptive and fixed binnings. 



F ioo 


Adaptive (r=2.0) 


Fixed (T=2.0) 


Adaptive (T=2.4) 


Fixed (T=2.4) 


5 


1.02+0.01 


1.03+0.01 


1.05+0.01 


1.03+0.01 


1 


1.00+0.02 


1.04+0.02 


1.05+0.02 


1.05+0.02 


0.5 


1.00+0.03 


1.03+0.02 


1.05+0.04 


1.01+0.05 



Table 5. Comparison of the average rms of the Af distributions for adaptive and fixed binnings. 



3.3.3. Interbin correlation 

In the adaptive-binning method, the starting time of a given bin depends on the starting time and flux-dependent width of the 
previous bin and thus by recurrence on the whole flux history since the start of the light curve. In this section, we explore the impact 
of the interbin correlation on the measured source parameters by comparing these in two consecutive bins. A comparison will again 
be made with the results of a similar analysis performed with a fixed binning, for which fluxes measured in consecutive bins are 
independent by nature. The flux (index) in one bin is plotted versus the flux (index) in the next bin in Figure [TT] left (right) for the 
whole light curve of a particular realization. 
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Fig. 9. Distributions of mean values (left) and rms (right) for flux (top) and photon index (bottom) obtained in adaptive- (solid) and 
fixed-binning (dashed) light curves. In this example, T = 2.4 and F\oo-5. 



The distribution of correlation factors obtained with the different realizations are presented in Figure [12] No significant differ- 
ences are observed between distributions obtained with adaptive and fixed binnings for all cases. Both distributions are centered 
close to (as expected for the fixed-binning case), demonstrating that the interbin dependence intrinsic to the adaptive binning 
method has little impact on the correlation between parameters measured in consecutive bins. 



4. High-level analysis with variable sources 

4.1. Duty cycle 

One way to describe variations in a light curve is via a duty cycle, corresponding to the fraction of time a source is in a "bright" state. 
More generally this can be presented as a flux distribution or its integral, the cumulative distribution function (CDF). In Figure[T3lwe 
have plotted flux, normalized to a mean of 1, versus 1-CDF, for the light curves displayed in Figure Q] The curves give the fraction 
of time (1-CDF) that a source is above a particular flux level. Four curves are shown for each light curve, the 'true' light curve, the 
adaptive-binning light curve, the fixed-binning light curve (with the same numbers of bins as the adaptive-binning light curve) and 
the fixed-binning light curve with higher (16 times) time resolution. The two fixed-binning light curves are complementary in that 
the high resolution one works well at the high-flux end but fails at the faint-end where the noise governs the distribution because 
of the low signal-to-noise ratio of the individual points. With adaptive binning the bin widths can be optimized such that in the 
presence of noise, the flux distribution, or duty cycle, can be described over a much wider range than is possible when a fixed bin 
width is used. The larger the range in fluxes the more advantageous the adaptive binning is for this type of analysis. 
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Fig. 10. Distributions of correlation factor between spectral index and flux for adaptive (solid) and fixed (dashed) binnings. T=2.0 
(left), T=2.4 (right), F m -5, 1, 0.5 from top to bottom. 
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Fig. 11. Fluxes (left) and photon spectral indices (right) measured in consecutive bins plotted one against another for a given 
realization. 



4.2. Power Density Spectrum (PDS) 

The PDS is one of the main tools used to characterize light curve variability. For evenly sampled data the PDS is computed directly 
via the Fourier transform of the light curve. For light curves with uneven sampling the Lomb-Scargle algorith m can be used and 
the PDS distortio ns caused by the sampling are usually modeled by comparison with simulated light curves dDone etal.lli99l 
Uttlev et al.ll2b 02). This method cannot easily be applied to adaptively binned light curves since the bin widths depend on flux, 
making time sampling different for each simulation. It is still possible to make a qualitative estimate of the PDS for an adaptive- 
binning light curve. To do this we oversample the light curve with a linear interpolation between the center of each of the adaptive 
time bins. The PDS is then computed with a standard Fourier transform. The binning and interpolation modify the power density 
around and above the frequencies corresponding to the widths of the time bins. To check at what frequencies this becomes important 
for a particular light curve we compute the distribution of Nyquist frequencies (equal to 0.5 /dt where dt is the bin width) associated 
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Fig. 12. Distributions of correlation factors between consecutive bins both for flux (top panels) and photon spectral index (bottom 
panels). r=2.0 (left-hand panels) and T-2A (right-hand panels), Fioo=5, 1, 0.5 from top to bottom in each panel. Solid: adaptive 
binning, dashed: fixed binning 



with the bin widths. In Figure [14] this distribution is plotted for the three simulated light curves shown in Figure [T] in each case 
together with the PDS for the adaptive- and fixed-binning light curves as well as for the true light curve. For the fixed-binning case 
we also show the PDS after subtraction of the white noise level, illustrating where, for this PDS, the signal is lost in the noise. It is 
clear from these examples that the PDS of the adaptive-binning light curve provides a qualitatively reasonable estimate of the PDS 
over a similar frequency range as the PDS for the fixed binning case. At the higher frequencies however, the adaptive-binning PDS, 
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1-CDF 1-CDF 

Fig. 13. The normalized flux is plotted against the parameter 1-CDF, obtained for the different light curves shown in Fig.Q] Left: 
whole range, right: zoom on the low-CDF region. Solid purple: true light curve, dashed green: fixed binning, dotted red: fixed 
binning with higher resolution, dot-dashed blue: adaptive binning. 

as we compute it here, is determined essentially by the binning interpolation. A more detailed discussion of the PDS properties and 
ways in which the binning effects can be taken into account is outside the scope of the present paper. 

5. Summary 

We have presented a method enabling constant-uncertainty/constant-significance light curves to be computed with the Fermi-LAT 
data. These light curves encapsulate more information than do fixed-binning light curves, without favoring any a priori arbitrary 
timescale and while avoiding upper limits. Although primarily developed for blazar studies, the method is applicable to all variable 
sources. The current implementation allows the time intervals to be calculated in a simple and quick way, while the reported 
parameters result from the standard LAT analysis. A slight skew in the measured flux distribution is a feature of the method. 
Possible adverse effects were explored by Monte-Carlo simulations but no significant issues have been revealed. This method opens 
capabilities regarding timing analysis in the y-ray domain that were only available so far at longer wavelengths. 

6. Acknowledgments 

It is a pleasure to thank S. Fegan and J. Scargle for useful discussions and suggestions. 

The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have sup- 
ported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics 



12 



B. Lott et al.: Adaptive-binning method for Fermi-LAT light curves 



and Space Administration and the Department of Energy in the United States, the Commissariat a FEnergie Atomique and the Centre 
National de la Recherche Scientifique / Institut National de Physique Nucleaire et de Physique des Particules in France, the Agenzia 
Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and 
Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) 
in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden. 
Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di 
Astrofisica in Italy and the Centre National d' Etudes Spatiales in France. 



13 



B. Lott et al.: Adaptive-binning method for Fermi-LAT light curves 



>r 2 



O 



-2 



CD 
O 



1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 


1 1 1 1 1 1 1 


a ' 








: □ ..niw. 


<!!! ! ff 
I'M 





-3.0 



a- 2 



o 



-2 



-3.0 



-2.5 



-2.0 -1.5 -1.0 -0.5 
LOG FREQUENCY (day 1 ) 



0.0 



0.5 



1 1 1 1 1 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 


b : 




V 

V/l l V' V -/ i ' , .t r 








•if : 




i 1 1 i «1 * 

' ' '! 1 A 


4 ; 


.... 1 




.IN .' 



-2.5 



-2.0 -1.5 -1.0 -0.5 
LOG FREQUENCY (day 1 ) 



0.0 



0.5 




-2.0 -1.5 -1.0 -0.5 
LOG FREQUENCY (day 1 ) 



0.5 



Fig. 14. Comparison of Power Density Spectra for the same simulations as in Fig.Q] In each plot the solid curve (purple) is the PDS 
of the true light curve. The dashed (green) and long dashed (red) curves are the PDS calculated for the fixed binned light curves, the 
white noise level being subtracted for the latter. The PDS of the adaptive binned light curve is plotted as a dash-dot curve. All PDS 
curves were smoothed by averaging over logarithmic frequency bins. The distribution of Nyquist frequencies corresponding to the 
time widths of the adaptive bins is shown by the histogram (blue). Note that this distribution is plotted in linear scale, the zero level 
corresponding to the lower limit of the vertical axis. 
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Appendix A: 

The goal is to estimate the relative integral-flux uncertainty cr ln F from the data given particular source and background models. 
Both flux and cr\ a p estimates are based on an unbinned maximum-likelihood method presented in the following. These estimates 
have been found to be in agreement with those one s given by the pylik elihood analysis within a few percent. Some of the formulae 
derived below can also be found in the appendix of Abdo et al. (2010al). 

A.1. Log Likelihood 

The likelihood X. is the probability of obtaining the data given an input model which is in this case the distribution of the y-ray 
fluxes over the sky. The maximum likelihood method con sists in ma ximizing £, over the parameter phase space of the model to get 
the best match between the model and the data. Following ICashl(ll979l) fEa. 7), the lo garithm of the likelihood is obtained by adding 
up the contribution of each photon, with index i in the following, within a region of interest (ROI) over a time interval [To, T\]: 

rn£=^lnM i -jV pred (A.l) 

i 

where M, = Exp(Ei,To,Ti)[S s (Ei)PS F(r,, Ej) + Sb(Ei)] is the model event density, £, is the photon energy, r, is the angular 
deviation from the source position, Exp(Ei, To, T\) is the exposure (time x effective area) over the considered time period, S s (Ej) is 
the differential source spectrum, assumed to be a power law in the following, 

Er 1 



S ,(£,-)= A \^-J (A.2) 

where Eo is a reference energy, A is the differential flux for Ej = Eo and F is the photon spectral index, PS Fiji, E,) is the point 
spread function^, S b(E,) is the differential background per solid angle unit. Finally, 



AVd = j Exp{E,To,T x ) S S (E) + j S B (E)d£l 



dE (A3) 



is the number of photons within the ROI predicted by the model. 
A.2. Significance/Flux calculation 

The difference between /«£ assuming that the source is present or absent is: 

A In £ = In [1 + gin, Ed] - N src (A A) 

i 

where g(r^ Ej) is the local source/background ratio 
S s (E,)PSF(r h E,) 

gin, Ed = — — (A.5) 

JB(Ej) 

and N src is the predicted number of photons coming from the source only: 

N SIC = J Exp(E,T ,T l )Ss(E)dE (A.6) 

The flux of the source is estimated by maximizing A In X with respect to A (T is usually set fixed in Step 1). 

The Test Statistic (TS) is defined as 2Aln£. In principle, TS behaves as ax 2 distribution with a number of degrees of freedom 
equal to the number of free source parameters, but in our particular case the statistic is modified as negative flux values are not 
allowed in the fit ("1/2 x 2 ")- For large enough one has: TS - N%.. The TS of the source is the sum of the contribution of each 
photon within the ROI with a weight depending on r,. As long as the ROI radius is large enough to include all photons from the 
source, the TS computed over the ROI is independent of the ROI size. Indeed, the contributions from distant photons are negligible 
since there are suppressed by the Point Spread Function. 

A.3. Relative flux uncertainty 

The relative flux uncertainty cr\ a f is calculated in Step 1 using the formulae below so as to consistently estimate cr\ n f as obtained 
in Step 2. These formulae are derived under the assumption that F is left free in the fit performed in Step 2 (otherwise, one would 
just have cr ln F = x)- 

The elements of the Hessian matrix, the inverse of the covariance matrix, can be calculated as: 
ff„ = dHn£ - d2Mn£ (A.7) 



assumed to be independent of the angle between the satellite axis and the arrival direction of the photons. This is justified for the LAT. 
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where the parameters x\£ are either A or T. 

The relative statistical uncertainty on the differential flux is: 



AA \Haa\~ 1 ' 2 



^ (1 + 



gin, Ed 2 



(1+ g(n, Ed) 2 



■1/2 



The cross-term element Hat is: 



1 

Hat — — r 



I 



g(r i ,£ i )ln(£ j /£ ) 



-JOT 



ExpiE, T Q , T l )SsiE)PSFir,E)lniE/E )d£ldE 



(l+gi^Ed) 2 

On average (over many observations of the same source), Hat can be written as 

,E) 2 



Exp(E, T , T,)S B iE)\n(E/E 



gir,E) 



-dQ. 



(A.8) 



(A.9) 



(A. 10) 



assuming that the model is correct, i.e., the actual photon energy distributions are given by S siE) and S siE) for the source and 
background respectively. Defining the pivot energy E p j VOt as the value of the reference energy Eo (Eq. IA.2l ) for which Hat = 0, one 
obtains: 



pivot) 

where 



f WiE) In EdE 
J WiE)dE 



W(E) = ExpiE,To,T x )S B (E) 



r,E) 2 



gir,E) 



-dQ. 



(A. 11) 



(A. 12) 



In the following, the reference energy Eo is set to E p j VOt , as AA/A is minimum in that case. The Hessian matrix is then diagonal and 
thus so is the covariance matrix. 

F is the integral flux between E m - m and +oo. It can be written as: 



-L 



Ss (E)dE = ^lp± 

E„,i„ 1—1 \fc pi vo t 



\-T 



(A. 13) 



The relative uncertainty of F, cr\ n F , is obtained by propagating the errors of A and F, making use of the fact that the covariance 
matrix is diagonal for this choice of Eq : 



0"ln F 



I + crj; | In 



E pivot 



l 



r- 1 



1/2 



(A. 14) 



where cr r is the uncertainty of T, <x r = \H yy \ l ^ 2 



A.4. Optimum energy 

Using the optimum energy E\, corresponding to the energy for which the second term in Eq. lA.14l cancels. for E m i„ will lead to the 
minimum cr ln F relative to other choices of E m - m \ 



In E x = In E pivo , - - — 7 (A. 15) 



1 

r~ l 

In that case <r ln F is equal to AA/A and is obtained from Eq. lA.8 
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